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I. INTRODUCTION 

Derivation of macroscopic transport laws ab initio 
from microscopic laws of motion is still one of the hot top- 
ics in modern mathematical physics. The main question 
usually investigated is what type of transport emerges 
from microscopic laws, the two extreme cases being 
ballistic or diffusive transport. Yet, such classification 
seems elusive, especially in the field of strongly correlated 
quantum systems. To nevertheless explain the emer- 
gence of macroscopic transport laws, phenomenological 
approaches, viable under certain assumptions, are often 
used. One such approach is based on the well-known 
kinetic theory of gases,t3 devised by Boltzmann, where 
macroscopic transport laws emerge due to (quasi)particle 
scattering. Also routinely used is linear-response theory, 
which gives a direct method of calculating transport coef- 
ficients using variants of Green-Kubo formulas^ assum- 
ing local quasiequilibrium. Our approach here is differ- 
ent. Starting from the equations of motion (Schrodinger 
equation or master equation) we use extensive numeri- 
cal simulations to study transport in the one-dimensional 
Heisenberg model. 

Low-dimensional spin models have been studied from 
the very beginning of quantum mechanics. The Heisen- 
berg model, suggested'^ by W. Heisenberg in 1928, was 
initially proposed to explain a high phase transition tem- 
perature in ferromagnets that could not be accounted for 
by any other known direct interaction. Interaction be- 
tween nearest neighbor atoms, described by the Heisen- 
berg model, the so-called exchange interaction, is an ef- 
fective one and comes about due to the Pauli exclusion 
principle. Another important source of motivation to 
study such simple one-dimensional (ID) models, in par- 
ticular antiferromagnetic ones, conies due to the fact 
that they represent simplest models of strongly corre- 
lated electronic systems, much studied in past decades. 
In addition, one-dimensional spin syst ems are realized in 
real so-called spin-chain materials.'^ Experiments have 
shown, see e.g. Ref. [SI that ID spin chains found in 



such materials, for instance the isotropic Heisenberg one, 
have a pronounced effect on transport properties, giving 
additional boost to theoretical studies. A great deal of 
research was devoted to the transport properties of the 
anisotropic Heisenberg model {XXZ model) . Although it 
is one of the simplest ID spin models, being even solvable 
by the Bethe ansatz,'^ there are still many open questions 
concerning its transport properties, including whether 
finite-temperature spin transport is ballistic or diffusive. 
Classification of transport regimes was the main motiva- 
tion for our work. 

Let us briefly review known facts about transport in 
the ID Heisenberg model; for more extensive reviews see 
Refs. [51IH1 andP, As forementioned, one of the standard 
approaches for studying transport properties of quantum 
systems is linear-response theory. Ballistic and diffusive 
transport can be distinguished via an observation of the 
Drude weight, the prefactor of a 5 function at zero fre- 
quency in the frequency-dependent transport coefficient. 
Nonzero Drude weight signals ballistic transport in the 
linear-response regime. The question of energy trans- 
port in the XXZ model is simple: energy current is a 
conserved quantitj'^^'^ and therefore energy transport 
is ballistic. For the dependence of the thermal Drude 
weight on parameters, see Ref. [TT] and references therein. 

In the present work, we shall focus on magnetization 
(spin) transport, which is much less understood, with 
only few rigorous results. It has been shown that the 
spin Drude weight is nonzero (i.e., magnetization trans- 
port is ballistic) at zero temperature^ in the gapless 
phase for A < 1 as well as at infinite temperatur^i^ 
for A < 1, where A is the anisotropy. It is reasonable 
to expect, a nd al so supported by quantum Monte Carlo 
calculations j ^ 1 that transport is ballistic also at finite 
temperatures. On general grounds, a lot of attention 
has been devoted to the connection between integrabil- 
ity and the nature of transporlP^ being either ballistic or 
diffusive. A recent solvable diffusive modeP^ shows that 
solvability does not necessarily imply ballistic transport. 

While the spin transport in the gapless phase is rel- 
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atively well understood, the behavior at the isotropic 
point A = 1 and in the gapped phase A > 1 is hotly 
debated. The main difhculty is that numerically it is 
very hard to access behavior in the thermodynamic limit, 
while there are only few analytical approaches. A no- 
table one is in terms of Mazur's inequalitjff^, which can 
be used to bound the Drude weight away from zero if 
a conserved quantity exists that has a nonzero overlap 
with the magnetization current.^ Unfortunately, in the 
half-filled case (zero total magnetization) and A > 1 
of interest here, no such quantity is known. Numeri- 
cal methods like exact diagonalizatiorfi^tsni or the Lanc- 
zos methocPi' are all limited to small systems of few 
10 spins, making thermodynamic extrapolation difficult. 
A relatively recent method is a time-dependent density- 
matrix renormalization-group (tDMRG) procedure that 
enables simulation of ID nearest-neighbor systems of sev- 
eral 100 spins. It has been used successfully to study the 
spreading of wavepackets in the XXZ model. Particu- 
larly useful is its master equation variant, where one has 
a genuine nonequilibrium setting, enabling one to study 
also far from equilibrium situations. It has been used 
to show a diffusive tran sport in the gapped phase at in- 
finite temperaturejI^^ESl which has been also confirmed 
using correlation function^^H and the projector operator 
method. Analytical studies of the master equation de- 
scribing nonequilibrium XXZ model in the gapped phase 
frequently have difficulties. One problem is that pertur- 
bative treatments often have zero convergence radius in 
the thermodynamic limit, for instance, a perturbative se- 
ries in the coupling to the reservoir^i^ for A > 1 or a per- 
turbative serie^^H /^^ qj. they have a finite convergence 
radius but are difficult to obtain in the thermodynamic 
limit, as is the case for large A in the gapped phase.l^ 
Perturbative studies in 1/A sug gest t hat the diffusion 
constant decreases as A increases.'^fiES 

Because most spin-chain materials realize the isotropic 
Heisenberg model, the point A = 1 is of particular inter- 
est. It is also the most controversial one. Analytical 
Bethe ansatz calculations give contradictory results, in- 
dicating zerd^ or a nonzero Drude weight,'^ with the 
problem being how to properly account for all states 
important at a nonzero temperature. Quantum Monte 
Carlo calc ulation s,^*^ ^'"^ bosonization,^'^ and exact diag- 
onalizatiorPsnEU predict a nonzero Drude weight at fi- 
nite temperatures for A = 1. Based on bosonization'^- 
and numerically calculated current autocorrelation func- 
tion using a tDMRG method, a zero (or small) Drude 
weight at nonzero temperature is advocated in Refs. EH 
and [37l whose results are also supported by quantum 
Monte Carlo calculation in Ref. HU] A recent result,'^ on 
the other hand predicts anomalous magnetization trans- 
port at infinite temperature and A = 1, with the diffusion 
constant diverging as ~ VL with the system size. 

An interesting future possibility to study ID strongly 
correlated systems is via controlled experiments with cold 
atoms. Experimental quantum optical techniques have 
advanced to the point where it is possible to realize such 



models in a controlled environment of optical latices or 
ion traps. An advantage of such an approach is that 
one can choose the values of system's parameters at will. 
First realizations of exchange interaction or of simple spin 
systems have already been achieved.'^Hl 



The main goal of the present paper is to study the 
magnetization transport at finite temperatures in the 
gapped phase of the anisotropic spin-1/2 Heisenberg 
model. Two numerical approaches, both based on the 
tDMRG method, will be used: one is based on observa- 
tion of time evolution of magnetization profiles for initial 
nonequilibrium pure states, while the second one is based 
on studying the magnetization current in the nonequilib- 
rium steady state of an open quantum system described 
by a quantum master equation. By the first method, 
we could, in principle, discriminate between the ballis- 
tic or diffusive behavior by comparing the evolution of 
expectation values of magnetization (z-spin component 
at each chain site) to that expected from macroscopic 
transport laws, provided we would be able to simulate 
very long chains for a very long time. Unfortunately, 
this is not the case and the results for pure state evolu- 
tion are rather inconclusive, showing a mixture of ballis- 
tic and diffusive characteristics. In the master equation 
approach though, performed at higher energy densities, 
one can give a quantitative prediction about the trans- 
port by studying the scaling of the magnetization cur- 
rent with the system size at a fixed driving. The two 
methods work best in the complementary temperature 
regimes. The one for pure states is best at low tem- 
peratures, where a state only locally deviates from the 
ground state and its entanglement is small, while the 
master equation simulation with density matrices works 
best at high temperatures where the operator-space en- 
tanglement of a density matrix is small. Both methods 
have been used before to study the magnetization trans- 
port in the XXZ model, pure-state method in Ref.lMland 
master equation in Ref. 1251 however not in the temper- 
ature regime considered in the present paper. We also 
point out that with a pure-state evolution at very low 
temperatures one is not able to access the asymptotic 
transport regime with present computers, so some care 
has to be taken making statements about the transport. 



The structure of the paper is as follows. In Sec. |TT] 
the Heisenberg XXZ model is defined and the tDMRG 
method for evolution of pure and mixed states is briefiy 
described. In Sec. |III[ main results concerning transport 
properties are presented, with the analysis of the evolu- 



tion of pure initial Gaussian-shaped states in Sec. IIIB 



and the results from the master equation setting in 
Sec. nil CI In Sec. IIVI the results of the evolution of do- 
main wall-states are presented. 
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II. MODEL AND METHODS 
A. Model 

A one-dimensional Heisenberg XXZ model is defined 
by the Hamiltonian 



H 



L-l 

1=1 
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(^+5r,i+H.c.)+A^f5, 



(1) 



with spin-1/2 operators Sf ^ Sf , S[ for l-th position on 
chain, and = Sf ± iSf is spin raising/lowering oper- 
ator. The spin chain is composed of L spins. Coupling 
strength will be fixed at J = 1, and open boundary con- 
ditions will be used. A is the only remaining parameter 
and determines the anisotropy of the model. 

The above model, written in the spin language, can 
be transformed to a model of spinless fermions using the 
Jordan- Wigner transformation,!^ resulting in a Hamilto- 
nian 



L-l 



1=1 



^( t 



H.c.) + A(nz-i)(n,+i-i) 



(2) 

where c;, c| are standard fermionic annihilation/creation 
operators, while n/ — cjci. Magnetization (spin) trans- 
port in a spin chain given by Eq. JTl) is thus equivalent 
to a transport of particles in Eq. (pE Expectation value 
of Si is directly linked to fermionic particle density as 
ni = Sf + 1/2. Both descriptions can therefore be used 

interchangeably. Total magnetization M = X^^Li i'^i) 
is conserved quantity in the XXZ model; in fermionic 
picture, it corresponds to the conservation of the total 
number of particles, n = X^z^i 



B. Methods 

We have used two closely related variants of the 
time-dependent density-matrix renormalization-group 
(tDMRG) method: one for the evolution of pure quan- 
tum states and another for the evolution of a density 
matrix describing a system coupled to reservoirs. Both 
methods are based on writing expansion coefficients of 
a state in terms of products of matrices, the so-called 
matrix product state (MPS) ansatz for pure states and 
matrix product operator (MPO) ansatz for density ma- 
trices. In the following, we will only briefly overview the 
tDMRG method. For a more complete presentation see 
original references.'^ 

For the determination of a profile evolution from a 
given initial state, a time-dependent Schrodinger equa- 
tion has to be solved. 
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State vector ji/'), spanned by the tensor product of local 
Hilbert spaces for each site, can be written as 



(4) 



where s enumerates all basis vectors of the Hilbert space 
s — (si,...,sl) with si G {0,1} for the spin-down/up 
states. The dimension of Hilbert space for the spin chain 
of length L is 2-^. Solving the Schrodinger equation using 
exact diagonalization is thus possible only for very small 
systems. However, it turns out that often not the whole 
Hilbert space is relevant for the solution and much larger 
systems can be solved by a clever choice of basis vectors, 
for instance, as is done in the tDMRG method. In MPS 
formulation, expansion coefficients Cg are expressed as 
traces of product of L matrices A^', I — 1,... ,L, of 
dimension K x K, 



= ti:{Al^---A2). 



(5) 



Time evolution of the system is therefore described by 
time-dependent matrices Aj*' (t) , which can be efficiently 
calculated if the propagator U{t) — exp {—iHt) can be 
factorized into a product of unitaries that act only on two 
nearest-neighbor sites of a chain. This can be approxi- 
mately done for short time-step propagation U{t) using 
Suzuki- Trotter expansion. The system Hamiltonian is 
separated into two parts, H = Hi + H2, where all terms 
grouped inside each Hi and H2 mutually commute. Gen- 
eral Suzuki- Trotter expansion can be written^^ as U{t) — 
Hfe exp(-iafciJiT) exp(-«/3fciJ2T) + 0{tP+^), where the 
order p depends on the actual scheme used. The 0(tP+^) 
contribution gives rise to Suzuki- Trotter error, which can 
be reduced by either using expansions of higher order 
(having larger p) or reducing the step size r. After each 
time step, the MPS form of the state has to be restored 
using the SVD algorithm. In the process the matrices 
are enlarged and to prevent an exponential growth 
of their dimension they must be truncated back to dimen- 
sion K. This produces a truncation error etrunc which de- 
pends on the entanglement of the state being described. 
A cumulative truncation error due to truncations at each 
time step can serve as a very rough estimate of the pre- 
cision. A more robust method to check accuracy though 
is to simply increase K and check that the results do 
not change. By using imaginary time step r — zr, the 
tDMRG method can also be used for obtaining a ground 
state of a given Hamiltonian. 

The tDMRG method for evolution of pure-states can 
be extended to the evolution of mixed states by intro- 
ducing a 4^-dimensional Hilbert space of operators, with 
arbitrary operator given by 



ip) = E 



Cs a- 



where a- 



(3) {0,1,2,3}, with (jO = 1 



,sl), 



and 
^3 . 
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while lower indices denote a site along a chain. Follow- 
ing the reasoning of the MPS formulation for pure state 
evolution, coefficients can be written in a MPO form, 
analogous to the MPS ansatz in Eq. ([s]), with s/ enumer- 
ating the basis of single-site density operators. The time 
evolution of density matrix p is governed by the Lind- 
blad equation, with a formal solution p(t) = exp(£t)p(0), 
where £ is a Liouvillian superoperator. Thus the propa- 
gator for time step r can be again factorized as a prod- 
uct of local propagators using Suzuki- Trotter expansion, 
providing an efficient numerical scheme for the time evo- 
lution of a density operator. Because C also contains 
dissipative terms due to a coupling with reservoirs, some 
care must be taken to ensure that one has a Schmidt- 
decomposed form of MPO at each step. Details of our 
implementation can be found in the appendix of Ref. 1471 
In the simulations of pure states, the Suzuki- Trotter 
expansion of second order was used (p = 2), mostly with 
the time step of size r — 0.05. The dimensions of the 
decomposition K were chosen such that the truncation 
error on each time step during the evolution did not ex- 
ceed etrunc = 10~^ (required dimensions of decomposition 
K were of the order of ~ 100). For master-equation sim- 
ulations a fourth-order method with a time step r = 0.05 
has been used. The results were also verified by repeat- 
ing simulations at somewhat higher decomposition sizes 
K and smaller time steps, showing no deviation in re- 
sults. 



III. TRANSPORT PROPERTIES 

The main goal is to study magnetization transport in 
the Heisenberg model at low energies, particularly for 
A > 1, where magnetization transport looks diffusive at 
an infinite temperature. ^^^^^ We shall use two methods, 
one will be spreading of localized packets and the other 
the dependence of the spin current on the system size in 
the stationary nonequilibrium state of a chain coupled to 
different reservoirs at chain ends. It is useful to have a 
"thermometer" with which we will be able to determine 
the temperature to which respective simulations corre- 
spond. In the thermodynamic limit of stationary states 
of a master equation under weak driving the state is in a 
local quasiequilibrium and the temperature is a well de- 
fined concept. In the wavepacket simulations though, in 
which packets are initially far from equilibrium, we pre- 
fer to speak about energy density instead, as there is no 
local quasi-equilibrium and the temperature is not well 
defined. Another important scale in the gapped regime 
of A > 1 is the size of the energy gap between the ground 
and the first excited state. 

In this section, we shall therefore first establish the val- 
ues of the ground-state energy, the gap, and the relation 
between temperature and the energy density for the XXZ 
Heisenberg model, in particular for A = 1.5 used later. 
Next, we shall present the results for the wave-packet 
simulations, followed by master-equation results. 



Ground state energy, the gap, and the 
temperature 
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FIG. 1. Dependence of the canonical energy density Et Q 
on the temperature for the anisotropic Heisenberg model, 
A = 1.5, and open boundary conditions. Curves are for 
L = 8, 16, 32, 64, bottom to top, while horizontal lines de- 
note the average energy densities of Gaussian packets we use 
to study transport (note that all are much higher than the 
gap). 

In order to be able to estimate where in the en- 
ergy spectrum our initial states are, we have calculated 
the ground-state energy of the anisotropic Heisenberg 
model,'^ the size of the gap, and for both quantities also 
finite-size corrections. Using numerically exact diagonal- 
ization for small sizes and an imaginary-time tDMRG 
for longer chains (up to L = 64), we have determined 
by fitting that for anisotropy A = 1.5 and open bound- 
ary conditions the energy density in the ground-state 
(which is always from the sector with M = 0) is 



ho = 



(H) 



L - 1 



-0.5234- 



0.29 

L-1' 



(7) 



Here (H) is the expectation value in the ground state 
of H. Finite-size correction to the asymptotic energy 
density is therefore of the order 0{1/L). At L = 200, 
used in our simulations, the ground-state energy den- 
sity is ho ~ —0.5247. Interesting to note is that 
for periodic boundary conditions finite size correction 
is smaller, namely, the ground-state energy density is 
h^^"^ « -0.5234-0.89/i2. The gap between the ground 
state and the first excited state (within the same symme- 
try class) is for A = 1.5 and open boundary conditions 



El 



Eo ~ 0.068+ -. 

1j 



(8) 



For a chain of length L = 200 the gap is Ei-Eo- 0.102. 
As we shall see, the energies of our initial conditions will 
always be significantly above the ground state gap. At 
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the energy scale below the gap, transport is trivially insu- 
lating and we are not interested in this regime. We have 
also calculated the relation between the thermodynamic 
temperature and the energy density Et in a canonical 
state. 



Er, 



L-1 '■ 



Pt 



exp i-H/T) 
trexp(-i?/T)' 



(9) 



Results are in Fig. [T] Such relation can be used as a 
"thermometer'^'^ for a given local energy density Ei^ 
one can determine to what thermodynamic temperature 
this corresponds by equating Ei — Ex and solving for 
the temperature T. One should be aware though that 
the validity of the canonical distribution is by no means 
granted for integrable systems. In fact, for specially cho- 
sen local reservoirs within the Lindblad master equation, 
deviations from the canonical distribution can be signifi- 
cant for integrable systemsJ^Hl for discussion of open sys- 
tems with general nonlocal coupling to reservoirs, see, 
e.g., Ref.lHl 



B. Localized packets 

One way to study transport is to initiate a spin chain 
in an initial state that is nonequilibrium with respect to 
the Hamitonian generating time evolution only in a small 
localized region. In other words, one prepares a localized 
initial packet and then studies how such a packet spreads 
in time. To obtain a state that is out of equilibrium 
only in a small region of space, we put the chain into a 
spatially inhomogeneous external magnetic field Bi and 
use the ground state of such a system as an initial state 
for our evolution. 



Preparation of initial states 



To prepare the initial state we take the Hamiltonian 



Ho = H + J2S[Bi 



(10) 



1=1 



where H is the XXZ model Hamiltonian from Eq. ([T]). 
For such Hamiltonian, the ground state |^o) was obtained 
using imaginary-time tDMRG method. At time t = 
the magnetic field is then removed and the initial state 
evolves according to H. The actual spatial dependence 
of Bi must be chosen in such way that time evolution of 
magnetization enables us to distinguish between ballistic 
and diffusive behavior. Its detailed form will be given 
further on a per-case basis. 

Fermionic nature of the antiferromagnetic ground state 
of the XXZ model leads to Friedel oscillations in mag- 
netization profile (S'f ) which get more pronounced when 
either effective interaction between fermions (A) or mag- 
netic disturbance for obtaining initial state gets larger. 



Friedel oscillations have a characteristic wavelength of 
27r/(2fci7'), where kp = Tr/2. As intensity of these oscilla- 
tions can blur the effective magnetization profile, simple 
averaging over two neighboring sites was used in the ma- 
jority of calculations, S[ = ((S'f.i) -I- (S'f ))/2, similar to 
Ref. |24l The same averaging was also used for the energy 
density profiles Ei = {{Ei_i) + {Ei))/2, where the energy 
density operator is Ei = ^(S'/"S',+i + H.c.) + AS[Sj;^^. 

Spatial dependence of the initial external magnetic 
field for the Gaussian-shaped initial profiles is given by 



_ [l-(I- + l)/2]^ 

Bi^Be -Bo, 



(11) 



where B determines the intensity of an initial distur- 
bance, <7b its width and Bq an overall offset of magnetic 
field that is used to adjust the total magnetization M of 
spin chain. 



2. Evolution of magnetization profiles 

Instead of focusing on the growth of the packet's vari- 
ance with time, as for instance in Ref. we focus on 
the evolution of the whole magnetization profile. The 
reason is that, in the variance one can get spurious effects 
as we shall discuss at the end of this subsection in llll B 4l 
In this section, we always take A = 1.5. We first pick a 
moderate B = 2, packet width as — 5 and a compen- 
sating magnetic field Bq = —0.18, resulting in an initial 
state with M = 0. In Fig. [2) we show a density plot 
showing local magnetization along the chain for times up 
to t — 50. At few time slices, we also show magnetization 
profiles. Two similar plots are also shown for the energy 
density. The energy E = (?/'o|i?|'0o) of the initial state 
Itpo) is about 1.8 above the energy of the ground state, 
which is much more than the value of the gap that is 
El ~ Eq w 0.1. The average energy density of our state 
is ft- = E/{L — 1) = —0.516. In equilibrium, this would 
correspond to the canonical expectation at the temper- 
ature T « 0.17; see also Fig. [T] Note that we do not 
make any claim that the state is in or is close to being 
in local equilibrium, in which case one could use a local 
temperature. From the profile plots we can see that the 
two bumps, moving away symmetrically from the origin, 
spread ballistically. In the density plots these are visible 
as ballistic jets. The speed of the wave front is almost the 
same for the energy and magnetization spreading. Based 
on these results one would be tempted to conclude that 
the transport is ballistic. However, one should be aware 
that the statement about transport is an asymptotic one, 
that is, for long times. In relatively small chains avail- 
able {L = 200) it could well happen that we have not 
yet reached this asymptotic regime. In fact, looking at 
the magnetization profiles one can see that there is some 
nontrivial dynamics going on in the region between the 
bumps, and that the bumps change shape and height 
with time. Also, one can argue on general ground that 
the short-time dynamics for generic models is always bal- 
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ior might be very large for the small-energy wavepackets 
shown in Fig. [2] 



FIG. 2. (Color online) Time evolution of the initial Gaus- 
sian packet obtained with B = 2, as = 5 and Bo ~ —0.18 
(M = 0) for a chain with L = 200 spins and A = 1.5. Sub- 
figures (a.i) and (a.ii) show the evolution of local magnetiza- 
tion, and subfigures (b.i) and (b.ii) the evolution of the en- 
ergy density. In the gray-coded density plots (a.ii) and (b.ii) 
horizontal dashed lines denote times t — 0, 10, 20, 30, 40, 50 
at which cross-sections are shown on the subfigures (a.i) and 
(b.i). Two slanted lines in the spin density plot (a.ii) indi- 
cate ballistic spreading of the magnetization with the speed 

V ~ 1.53. Two dotted slanted lines in energy density plot 
(b.ii) indicate ballistic energy spreading with ve ~ 1.55. Both 

V and Ve were determined by fitting the peaks of spin and en- 
ergy profiles at various times. Two dotted horizontal lines in 
the plot of energy density profiles (b.ii) are the average en- 
ergy density h ~ —0.516 of the initial state and of the ground 
state ho ~ —0.525. In subfigure (a.i) the dotted line, over- 
lying t = 50 cross-section, indicates spin profile calculated 
with MPS decomposition size K — 125, demonstrating neg- 
ligible deviations from the spin profile calculated at A" = 90 
(underlying bold line) even at longest simulation times. 



listic. For an explicitly solvable quantum model, where a 
transition from a ballistic spreading at short times to a 
diffusive at large times can be shown, see Ref.|48l We will 
in fact see that master-equation simulations, presented 
later, indeed support such a scenario. They also indicate 
that the transition time to asymptotic diffusive behav- 
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FIG. 3. (Color online) Time evolution of the initial Gaus- 
sian packet obtained with B = 5 and as = 10, L = 200, 
A — 1.5. Subfigures (a.i) and (a.ii) show the evolution of lo- 
cal magnetization, and subfigures (b.i) and (b.ii) the evolution 
of the energy density. In density plots (a.ii) and (b.ii) horizon- 
tal dashed lines denote times t = 0, 10, 20, 30, 40, 50 at which 
cross-sections are shown on the subfigures (a.i) and (b.i). Two 
slanted lines in the subfigure (a.ii) indicate ballistic spread- 
ing of the magnetization with the speed v ~ 0.77, determined 
by fitting to the wavefront reaching Si — 0.05. Two dotted 
slanted lines in subfigure (b.ii) denote speed of energy spread- 
ing Ve ~ 1.55 as determined in Fig. [2] while solid slanted lines 
indicate speed of spin profile wavefronts from subfigure (a.ii) . 
Two horizontal lines in the plot of energy density profiles are 
the average energy density of the initial state h ~ —0.406 and 
of the ground state ho ~ —0.525. 



To nevertheless be able to better assess the nature 
of transport also in packet-spreading simulations, we go 
to higher energies. There is strong numerical support 
that th e spin transport at an infinite temperature is 
diffusive .'2SH23 xhe ballistic spreading of jets at low en- 
ergies, seen in Fig. [2j should therefore weaken, for in- 
stance, slow down, or even entirely disappear, if the 
transition time to diffusive behavior gets short enough 



at higher energies. To verify this hypothesis we have 
simulated packets at higher energies by simply increas- 
ing the amplitude B and the width aB of the initial 
magnetic field. Note that in this way we also produce 
states that are more strongly nonequilibrium, at least 
for short times. In Fig. [3] we used B — 5 and cr = 10 
{Bq = —0.86) to prepare the initial state in a sector with 
zero total magnetization. The energy density of such 
an initial state is h — —0.406 and is therefore signifi- 
cantly above the ground state (w 240 times the gap). 
This energy density would in equilibrium correspond to 
temperature T « 0.58. We can see from the results in 
Fig.[3]that the wavepacket front still spreads ballistically. 
There is one important difference though compared to 
the smaller-energy packet from Fig. [2} the speed of the 
magnetization front is smaller, while the speed of the en- 
ergy front is unchanged. Taking an even higher energy 
packet, obtained by B = 5 and as = 15 {Bq = —1.21), 
the speed of the magnetization front decreases even more. 
Results for such a packet having an average energy den- 
sity h « —0.352 (which would correspond to temperature 
T 0.74) are in Fig.jlj From the spreading of packets at 
higher energies, several things can be concluded. First, in 
the regime of times and distances studied there is still a 
ballistic front in the magnetization that spreads at a con- 
stant speed. This speed decreases as one increases the en- 
ergy, as it should in order to accommodate for a diffusive 
behavior at an infinite temperature. What happens with 
the front at larger times, for which one would need larger 
systems, is hard to infer. For the packet with the energy 
density h = —0.516 (Fig. [2]), the speed is w w 1.53, for 
the packet with the average energy density h w —0.406 
(Fig. |3]), the speed is w « 0.77, and for the packet with 
h « -0.352 (Fig. |4]), the speed is v a 0.55. The sec- 
ond observation is that the speed of energy spreading 
is different than the speed of magnetization spreading. 
The speed of the energy front is always approximately 
ve ~ 1.55, which is incidentally also equal to the speed 
of magnetization spreading at low energieJ^. Note that 
the energy transport is ballistic in the Heisenberg model 
because the energy current is a conserved quantity. The 
third observation, especially visible in the spreading of 
the energy density in Fig. |4j is that there is an energy 
dispersion. The shape of the energy packet changes with 
time. Parts with higher energy are "overtaken" by lower- 
energy excitations, visible as the stretching of the energy- 
density profile in front of the main peak. 

All these simulations of the spreading of localized pack- 
ets of magnetization show signs of ballistic spreading, ei- 
ther in terms of ballistic jets, or ballistically propagating 
wavefronts. One should be careful though about conclud- 
ing that the transport is also ballistic at large times and 
in the thermodynamic limit, as these features could ex- 
hibit a nontrivial long-time behavior. As mentioned be- 
fore, there could be a large time scale Tdiff , so that only for 
times larger than T^is transport starts to show purely dif- 
fusive character. One can in fact connect this time scale 
with the diffusion constant D. Assuming for instance an 




FIG. 4. (Color online) Time evolution of the initial Gaus- 
sian packet obtained with B = 5 and as — 15, L = 200, 
A — 1.5. Subfigures (a.i) and (a.ii) show the evolution of lo- 
cal magnetization, and subfigures (b.i) and (b.ii) the evolution 
of the energy density. In density plots (a.ii) and (b.ii) horizon- 
tal dashed lines denote times t = 0, 10, 20, 30, 40, 50 at which 
cross-sections are shown on the subfigures (a.i) and (b.i). Two 
slanted lines in the subfigure (a.ii) indicate ballistic spread- 
ing of the magnetization with the speed i; ~ 0.55, determined 
by fitting to the wavefront reaching Si = 0.05. Two dotted 
slanted lines in subfigure (b.ii) denote speed of energy spread- 
ing Ve ~ 1.55 as determined in Fig. [2] while solid slanted lines 
indicate speed of spin profile wavefronts from subfigure (a.ii). 
Two horizontal lines in the plot of energy density profiles are 
the average energy density of the initial state h ~ —0.352 and 
of the ground state ho —0.525. 



exponential decay of the time-dependent spin current au- 
tocorrelation function, and taking into account that D is 
equal to the integral of the autocorrelation function, one 
sees that the autocorrelation function decay time scales 
as oc D. Therefore, if one has a diffusive behavior with a 
diffusion constant D, this introduces a characteristic dif- 
fusive time scaling as Tdiff ~ -D. If diffusion constant D 
is very large (we shall see that this is likely the case), one 
can have different behavior than the asymptotic diffusive 
one for t <^ Tdiff. Our simulations of wavepacket spread- 



ing therefore cannot distinguish truly ballistic transport 
from a diffusive with a large diffusive constant. To distin- 
guish the two, one would have to look at the spreading of 
very wide and very shallow packets at long times, so that 
local deviations from equilibrium are small. Because with 
tDMRG one is limited to chains of few 100 spins this limit 
can not be attained with the present computational re- 
sources. To nevertheless be able to say something about 
the spin transport in the anisotropic Heisenberg model 
at low energies, we also performed tDMRG simulations 
of transport in the master equation setting, where we 
can simulate systems at higher energies and some of the 
above mentioned problems do not appear. 

Before doing that, we shall briefly present two interest- 
ing observations about the spreading of localized packets 
that are not directly related to transport. Readers inter- 
ested only in transport properties can skip the next two 
subsections and jump directly to Sec. [Ill C[ 



3. Short-time wavepacket pinch 

As an interesting side-remark, not directly related to 
our study of spin transport, at short times two wider 
packets in Figs. [3] and |4] undergo a pinch. Their central 
part first contracts at very short times, and only then 
begins to spread. This can be seen in more detail in 
Fig. [5] which shows short-time behavior of magnetization 
for the packet from Fig. |4] obtained with B = 5 and 
ctb = 15. For times smaller than w 20 the central part 
of the packet shrinks while the shoulders expand. 




FIG. 5. (Color online) Short-time pinch of the initial packet 
obtained with B — 5 and as = 15 (same data as in Fig. |4|. 
After releasing the external magnetic field, the top of the 
packet first undergoes a contraction before it begins to spread 
at later times. 



4- Magnetization offset Bq 

In this subsection we would like to make few comments 
about the influence of total magnetization of the spin 




FIG. 6. (Color online) Time evolution of the initial Gaus- 
sian packet obtained with B — 2 and as = 5 without a 
compensating homogeneous magnetic field, Bq — 0, thus hav- 
ing nonzero z-spin component M = 2. Subfigures (a.i) and 
(a.ii) show the evolution of local magnetization, and subfig- 
ures (b.i) and (b.ii) the evolution of the energy density. In 
density plots (a.ii) and (b.ii) horizontal dashed lines denote 
times t = 0, 10, 20, 30, 40, 50 at which cross sections are shown 
on the subfigures (a.i) and (b.i). Two slanted lines in the spin- 
density plot (a.ii) indicate ballistic spreading of the magne- 
tization with the speed v « 1.55. Two dotted slanted lines 
in energy-density plot (b.ii) indicate ballistic energy spreading 
with ve ~ 1.55. Both v and ve were determined by fitting the 
peaks of spin and energy profiles at various times. Two dot- 
ted horizontal lines in the plot of energy-density profiles (b.i) 
are the average energy density of the initial state h « —0.507 
and of the ground state ho ~ —0.525. 



chain M on the evolution of spin profiles, specially on 
the variance (T^(t) of the profile. As M is a conserved 
quantity of the XXZ model, it is determined by the ini- 
tial state, i.e., the values of the magnetic field B and 
B() in Eq. (11). Remember that Bq is used to tune the 



value of the total magnetization while retaining approx- 
imately Gaussian shape of the initial magnetization pro- 
file. In the linear response, transport properties of the 
XXZ model are known to be related to the value of M. 



9 





50 






40 




b 


30 




b 


20 






10 












- ■ M = 

M = 2 

■■■■ at^ 



10 20 30 40 

t 

FIG. 7. Time dependence of spin profile variance a^{t) — ctq 
for initial states obtained with B — 2, as = 2 and total 
magnetization M = 0, and M = 2 (time evolution of cor- 
responding profiles is shown on Figs. |2] and [6|. Dotted line 
represents function at^ fitted to the variance of A4" = 2 case, 
demonstrating ballistic growth. The M = Q case on the other 
hand does not have simple time dependence of variance on 
attainable time scales. 

For states of nonzero magnetization, M ^ 0, ballistic be- 
havior was predicted on the basis of Mazur's inequahty as 
the spin current operator has a nonzero overlap with the 
conserved quantities of XXZ Hamiltonian.'i^ For M = 
this overlap is zero and no conclusion can be made, thus 
leaving the possibility of a diffusive transport. We gen- 
erated initial states by two alternative choices of initial 
magnetic-field parameters B and Bq. In the first case, to- 
tal magnetization of the initial states was tuned to M = 
by using an appropriate compensating magnetic field Bq. 
In the second case, initial states were generated without 
a compensating magnetic field, Bq = 0, thus resulting in 
a nonzero (but small) total magnetization M ^ 0. As 
system size approaches thermodynamical limit i — > oo, 
there is no difference between the initial states generated 
by the above alternative choices of Bq because Bq — >■ 0. 
Yet, for the attainable system sizes L ^ 200 used in 
the numerical simulations, the details of the initial state 
preparation can have an observable effect on the dynam- 
ics. In the following we define the variance of the spin 
profile (T^(t) and summarize our observations about the 
influence of the initial state preparation on the variance 
and overall evolution. 

Time dependence of variance, defined in the fermionic 
picture of the XXZ model as 

^'^^E(^-^)'-("'W)' ^^ = li2^{Mt)), (12) 
1=1 1=1 

was extensively studied in Ref. ,24, There, a ballistic 
spreading was established for A < 1, while for A > 1.5 
and M = a diffusive transport was advocated based on 
linear growth of at intermediate times. For M ^ 
it was showrP3 that one gets t^. Linear growth of 

variance in the case of diffusive transport follows from the 



macroscopic diffusion law for magnetization transport, 
dS''{l)/dt = -Dd'^S^il)/df. We were able to confirm 
the ballistic growth of variance for A < 1. For A > 1, the 
situation is more complex. Depending on the details of 
the preparation of the initial state, more complex dynam- 
ics emerges, which can result in a linear variance growth 
for intermediate timeJ^Sl (see Fig. [t] at t ss 10 — 25). The 
actual shape of the central magnetization packet though 
does not follow the behavior that is expected from the 
macroscopic diffusion equation as can be seen in Fig. [2j 
While in the case of A < 1 the details of the initial 
state preparation do not have pronounced effects on the 
evolution (either by observing (t) or the profile directly, 
both show a ballistic behavior), the effects are stronger 
in the case of A > 1. For states with Bq = (Fig. [6]) 
the profile evolution as well as the time dependence of 
the variance (Fig. [?]) seem purely ballistic, at least at at- 
tainable time scales. For the M = states on the other 
hand, while the profile evolution still contains ballistic 
features (Fig. [2| , the variance does not grow as ct^ cx t'^ 
but instead has a nontrivial time dependence. While this 
nonballistic growth of variance may be taken as a hint of 
transition to a diffusive behavior, a detailed analysis of 
spin profiles reveals that observed behavior can be at- 
tributed to the negative magnetization "bumps" on the 
outer edges of spreading disturbance as seen on Fig. [2] for 
times t > 30. As the value of the variance CT^(t) strongly 
depends on the profile values further away from the mid- 
dle of the chain, such bumps have a strong impact on the 
variance. The occurrence of the bumps in the profile can 
be attributed to the emergence of a nonequilibrium mag- 
netization plateau in the central region of the chain, be- 
tween the packets traveling in opposite directions. Due to 
the conservation of total magnetization M, this nonequi- 
librium plateau of magnetization is compensated by an 
opposite deviation of magnetization on the edge of dis- 
turbance. Variance of the profile therefore should not be 
taken as a sole criteria for the determination of transport 
regime. Thus we have focused mainly on the speed and 
the behavior of spin-profile wave fronts, emerging from 
the central peak, as these appear to be less dependent 
on the details of the initial state preparation. Note that 
we avoid making any definitive statement about the na- 
ture of transport from the wave-packet evolution only. 
These simulations will serve us only as a guide to master 
equation simulations that we present next. 

C. Master equation setting 

To induce a nonequilibrium steady state (NESS) we 
couple our spin chain to reservoirs. The coupling is de- 
scribed in an effective way via a set of Lindblad opera- 
tors acting on the first two and the last two spins. The 
Lindblad master equation describing the evolution of the 
density matrix i^iZffil 

^o^i[p,H]+£^'^{p)^C{p), (13) 
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where the dissipative hnear operator C is expressed in 
terms of Lindblad operators Lk , 

C^'^p) = {iLkP,Ll] + [L,,pLl]) . (14) 

To induce a NESS at a finite temperature we use Lind- 
blad operators that couple to the first and last two 
spins of the chain ~ the so-called two-spin bath. De- 
tails of the two-spin bath implementation can be found 
in Refs. |25] and |47l There are therefore 16 Lindblad 
operators at each end. They are chosen in such 
a way that they would induce a grandcanonical state 
~ exp {—H/Tl^r + Pl,rM) on these two spins in the ab- 
sence of Hamiltonian evolution by H. Because in our case 
the evolution by H is present it introduces some bound- 
ary resistance effects that also affect the efficiency of the 
method. Due to these boundary effects it is rather dif- 
ficult to cool the chain to very low temperature^^. We 
use reservoirs with the same temperature at the left 
and Tr on the right end of the chain, while the chemical 
potential is fiL — 0.1 at the left and pR — —0.1 at the 
right end. Because of this symmetric driving the average 
magnetization in the NESS is zero, as is also the en- 
ergy current. To calculate p{t), and therefore also NESS 
given by limt_>oo p{t), we use the tDMRG method with 
a matrix product operator ansatz. After long time p{t) 
converges to a stationary nonequilibrium state whose ex- 
pectation values then give us transport properties. Due 
to boundary resistance the temperature in the bulk of the 
chain is not the same as the imposed temperature of the 
"reservoir" Lindblad operators. To determine the actual 
temperature in the system in the nonequilibrium steady 
state we use the expectation value of the energy den- 
sity as a "thermometer" ^ equating it to the canonical 
one, and thereby determining the effective temperature, 
as described at the beginning of this section. 

Because one has to simulate evolution of density oper- 
ators instead of pure states, open system formulation is 
computationally more demanding than pure-state sim- 
ulation. This typically means that somewhat smaller 
chains can be simulated. In addition, simulation at low 
energy (temperature) is more demanding because the 
operator-space entanglement of the NESS increases with 
decreasing temperature.!^ For instance, at an infinite 
temperature the NESS, being proportional to 1, is sepa- 
rable with no entanglement. Therefore, due to computa- 
tional constraints, we had to focus on somewhat higher 
energies than in the wavepacket simulations. 

Once we determine NESS for a given driving and 
length L we calculate the expectation value of local 
magnetization, obtaining the difference in magnetiza- 
tion between left and right ends AS*^, local spin cur- 
rent ji = (S'fS'f,^;^ — SfS[_^^) (which is independent of 
Z), and local energy density. Because the deviation from 
the equilibrium zero magnetization is small, around 0.1 
in high energy simulations and « 0.01 at the lowest en- 
ergy, and the imposed temperature is the same at both 
ends, the energy density in the NESS is constant along 
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FIG. 8. Dependence of the diffusion constant obtained from 
NESSs on the inverse temperature. Dashed line suggests an 
exponential dependence on the inverse temperature, D « 
0.5 exp (1.6/3). On the top axis we also list temperatures and 
energy densities corresponding to 5 NESSs. Inset: scaling of 
finite-size D with the chain length L. Horizontal lines indicate 
the asymptotic D's used in the main plot. 



the chain. The effective temperature in the bulk of the 
chain is then determined by equating this energy den- 
sity to the canonical one (local chemical potential in the 
NESS is close to zero). The finite-size diffusion constant 
can then be determined a.s D = (L — 1) • ji/AS^. Taking 
into account definition of spin conductivity and of D cal- 
culated herejISSl the spin conductivity k can be obtained 
from the diffusion constant simply as k = /3 D {/3 — 1/T). 
To ensure that the behavior is really diffusive and such 
D does not depend on size L, we have calculated D for 
sizes up to L = 64. This data can be seen in the inset of 
Fig. [sj One can see that at high temperatures (T > 2, 
lower three chain lines in the inset) a good convergence of 
D is obtained, signaling diffusive spin transport also at a 
finite temperature. At the lowest temperature T « 0.82 
that we were able to simulate, convergence is less clear. 
From the data on D{T) in Fig.js] one can clearly see that 
the diffusion constant increases at lower energies, the in- 
crease being perhaps exponential in the inverse temper- 
ature /?, as indicated by a dashed line. At low tempera- 
tures/energies the diffusion constant can get very large. 
If we dare to extrapolate this dependence to the energy 
of the packet simulated in Fig. |2j and we take the average 
temperature T w 0.17 as a crude estimate, the diffusion 
constant would he D ^ 5000. This means that to really 
observe a diffusive behavior in a wavepacket simulation 
one would have to simulate chain up to times of the or- 
der ~ 5000, demanding also chains of a similar size. We 
expect that ballistic jet, visible at short times, will dis- 
appear after this long time scale. All these results mean 
that with the tDMRG, being limited to few 100 spins, 
one probably cannot conclusively say whether the spin 
transport is diffusive or ballistic at such low energies. 
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With an open-system tDMRG version we have nev- 
ertheless obtained strong indications that at not too 
low temperatures, T > 1, the anisotropic Heisenberg 
model is diffusive. The diffusion constant increases fast 
with decreasing temperature. Another thing we know is 
that at energy scales below the gap {Ei — i?o ~ 0.1 at 
A = 1.5), the system is trivially insulating; therefore, 
D{T = 0) = 0. At low temperatures that are still much 
above the gap there are then basically two possibilities: 
either the system is ballistic, meaning that the diffusion 
constant diverges at a finite temperature T, or the sys- 
tem is diffusive but with an exponentially large D. We 
find the latter scenario more plausible. 

In any case, as one decreases A toward 1, the gap dis- 
appears and therefore also an insulating state at T = 0. 
At A = 1 the diffusion constant is therefore infinite at 
zero temperature (in agreement with a nonzero Drude 
weighlP^, the transport therefore being either ballistic 
or anomalous. Recent results in Ref. [211 show that at 
an infinite temperature and A = 1 it is anomalous (su- 
perdiffusive) . 



IV. DOMAIN WALL DYNAMICS 

In the present section we will focus on particular initial 
states whose time evolution is quite different from the one 
of the Gaussian packets. These are states with a domain- 
wall-shaped initial magnetization profiles. Such states 
are nongeneric and therefore do not influence our con- 
clusions about magnetization transport reached earlier. 
The purpose is just to point out that, due to symmetry, 
there are particular states with different behavior. Such 
states we re inv estigated in the context of domain-wall 
dynamics^^'^ relaxation dynamics and domain- wall 
stability.'S3H56] Pqj. gapless regime (A < 1), numerical 
and analytical results suggest ballistic spreading of the 
initial domain-wall. Less is known about dynamics in a 
massive phase (A > 1). Time evolution of completely 
polarized domain walls have been studied numerically in 
Rcf. 63, and analytically using semiclassical approxima- 
tion in Ref. [6TJ Here we extend analysis to low-energy 
partially polarized domain states. 

We get domain-wall initial states as ground states of 
the Heiseneberg XXZ chain defined by Eqs. ([I]) and ( 10 ) 
in a step-shaped magnetic field given by 



Bi = 2Be{ 



L + 1 



l)-B, 



(15) 



where O is Heaviside step function and B is the mag- 
nitude of an external magnetic field that polarizes spins 
on the left and right side of the chain in opposite direc- 
tions. In the limit i? — > oo, ground state of the chain is 
given by the simple product state |t •.• ti •■• i); while it 
is more complicated for finite B. Nonetheless, the shape 
of the spin profile (S'f ) retains approximately step-like 
shape if we employ the nearest-neighbor averaging of S'f 
to account for Friedel oscillations. The main question we 




FIG. 9. (Color online) Evolution of spin profile of a partially 
polarized domain-wall-like initial state in a gapless regime, 
A = 0.5 and B = 0.5. In density plot (b), horizontal dashed 
lines denote times t = 0, 10, 20, 30, 40, 50 at which cross- 
sections are shown on the subfigure (a). Density plot (b) 
shows absolute value of spin profile \Si\. Ballistic spreading 
of the domain wall is clearly seen. 



studied was whether the initial domain-wall stays local- 
ized or it decays with time. To access this we have deter- 
mined the domain-wall width w at time t as the distance 
between locations, positioned symmetrically around the 
middle of the chain, where the value of magnetization S[ 
exceeds positive/negative offset value of magnetization 
±Mi„ equal to half of the maximal one. 

For the gapless case A < 1, the ballistic spreading 
of the domain wall was obtained by directly observing 
linearly increasing width of the domain wall {w ^ t). 
Ballistic behavior of domain- wall states was observed in- 
dependently of the level of polarization, determined by 
initial state polarization B. An example of profiles for 
B = 0.5, A = 0.5 is shown in Fig. |9] 

In the gapped regime A > 1, the dynamics of the ini- 
tial domain wall is more complex. To determine whether 
the domain wall spreads diffusively or remains frozen due 
to localization, we have observed time evolution of the 
domain- wall width w{t) for different values of A > 1 and 
for various polarizations of initial states, determined by 
B. For all A, the initial state used has been the same 
and was obtained as a ground state of Hamiltonian with 
A = 1. At t = magnetic field B was turned off and A 
was quenched to the target value. We have observed that 
for all A > 1, independently of the initial state polariza- 
tion B, the width of the domain wall w{t) grows until 
it attains some maximal value Wmax and then starts os- 



cillating around some final value iWfinai (see Fig. 10 for 
an example oi B — 0.5). Whether these oscillations are 
damped out as the time progresses could not be deter- 
mined on the accessible time scales. It therefore appears 
that for A > 1 the initial domain wall stays localized. 
The functional dependence of the domain-wall width on 
the value of anisotropy A was suggested in Refs. [55] and 
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FIG. 10. (Color online) Evolution of spin profile of a partially 
polarized domain-wall-like initial state in a gapped regime, 
A = 1.5 and B = 0.5. In density plot (b) horizontal dashed 
lines denote times t = 0, 10, 20, 30, 40, 50 at which cross sec- 
tions are shown on the subfigure (a). Density plot (b) shows 
absolute value of spin profile \Si\. Here the spreading of 
the domain wall stops and the domain wall starts to oscil- 
late around a stable shape (see, e.g., cross sections at times 
t = 40 and i = 50). 




FIG. 11. (Color online) Dependence of the maximal domain- 
wall width Wniax on the anisotropy parameter A for A > 1 
and the initial state of a completely polarized domain wall. 
Solid line has the functional dependence of Eq. ( 16 1 with the 
best fitting A = 1.12. 



states in the ferromagnetic Heisenberg system with kink 
boundary terms in Ref. 1651 Recently, localization of a 
fully polarized domain wall in the quantum anisotropic 
Heisenberg spin chain (ferromagnetic and antiferromag- 
netic) has been observed--- in the context of negative dif- 
ferential conductivity and explainecP^ in terms of a one- 
magnon localization. See also recent work in Ref. ISSl We 
have compared the measured maximal width^SH! w^^x to 
the suggested scaling form, Eq. (16). Results reported 
in Fig. [TT] are found to be in good agreement for large 
initial domain- wall polarizations {B — >■ oo), especially 
around Awl, where domain- wall widths span multi- 
ple spin sites. At larger anisotropics the domain wall 
widths are in a range of a few spins so Wmax is strongly 
dependent on the details of the procedure for the do- 
main width determination (e.g. averaging of magnetiza- 
tion profile and choice of domain- wall boundaries Af^,), 
so that the measured value of the domain width deviates 
from the suggested scaling form. For weaker polariza- 
tions of initial states {B ^ 1) the determination of the 
domain- wall width Wmax is more involved as transient ef- 
fects of the domain- wall dynamics get more pronounced. 
The functional dependence of Wmax for B 1 was there- 
fore not compared to the scaling form of Eq. ( 16 1. How- 
ever, we have noticed that the domain wall stays localized 
even for weakly polarized initial states, at least up to the 
timescales of t ~ 100. 



V. CONCLUSION 

We have studied magnetization transport at finite tem- 
peratures in the one-dimensional anisotropic Heisenberg 
model. By a combination of wavepacket spreading and 
the study of nonequilibrium steady states we reach a 
conclusion that the transport is diffusive in the gapped 
regime. By lowering the temperature the diffusion con- 
stant increases, perhaps exponentially so with the inverse 
temperature. A very large diffusion constant at low tem- 
peratures introduces a very long time and space scale 
that governs the transition from a ballistic behavior at 
short time to an asymptotic diffusive at long times. Using 
existing numerical techniques, that are limited to short 
times and small system sizes, it is therefore very difficult 
to observe diffusive behavior at very low energies. 



[5B^ for the ferromagnetic ID Heisenberg model, having 
the form 

w= (16) 

log(A + yA^^)' 

where A is dependent on the used definition of the 
domain-wall width. See also related interface ground 
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